Fit main model (see exploratory scripts and model comparison), visualize results.
library(tidyverse)
#> Warning: package 'tidyr' was built under R version 4.0.5
library(tidylog)
library(viridis)
library(marmap)
library(rnaturalearth)
library(rnaturalearthdata)
library(sf)
#> Warning: package 'sf' was built under R version 4.0.5
library(RColorBrewer)
#> Warning: package 'RColorBrewer' was built under R version 4.0.5
library(gganimate)
library(gifski)
library(latex2exp)
library(patchwork)
library(png)
library(RCurl)
library(wesanderson)
library(qwraps2)
library(ggcorrplot)
library(ggridges)
library(brms)
#> Warning: package 'Rcpp' was built under R version 4.0.5
# remotes::install_github("pbs-assess/sdmTMB")
library(sdmTMB) # sdmTMB_0.1.4.9004
# To load entire cache in interactive r session, do:
# qwraps2::lazyload_cache_dir(path = "R/analysis/condition_model_cf_cache/html")
# Specify map ranges
ymin = 52; ymax = 60; xmin = 10; xmax = 24
map_data <- rnaturalearth::ne_countries(
scale = "medium",
returnclass = "sf", continent = "europe")
# Crop the polygon for plotting and efficiency:
# st_bbox(map_data) # find the rough coordinates
sf::sf_use_s2(FALSE)
# https://stackoverflow.com/questions/68478179/how-to-resolve-spherical-geometry-failures-when-joining-spatial-data
swe_coast <- suppressWarnings(suppressMessages(
st_crop(map_data,
c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax))))
# Transform our map into UTM 33 coordinates, which is the equal-area projection we fit in:
utm_zone33 <- 32633
swe_coast_proj <- sf::st_transform(swe_coast, crs = utm_zone33)
# Define plotting theme for main plot
theme_plot <- function(base_size = 11, base_family = "") {
theme_light(base_size = base_size, base_family = "") +
theme(
legend.position = "bottom",
legend.key.height = unit(0.2, "cm"),
legend.margin = margin(0, 0, 0, 0),
legend.box.margin = margin(-5, -5, -5, -5),
strip.text = element_text(size = 9, colour = 'gray10', margin = margin(b = 1, t = 1)),
strip.background = element_rect(fill = "grey95")
)
}
# Define plotting theme for facet_wrap map with years
theme_facet_map <- function(base_size = 11, base_family = "") {
theme_light(base_size = base_size, base_family = "") +
theme(
axis.text.x = element_text(angle = 90),
axis.text = element_text(size = 7),
strip.text = element_text(size = 8, colour = 'gray10', margin = margin(b = 1, t = 1)),
strip.background = element_rect(fill = "gray95"),
legend.direction = "horizontal",
legend.margin = margin(1, 1, 1, 1),
legend.box.margin = margin(0, 0, 0, 0),
legend.key.height = unit(0.4, "line"),
legend.key.width = unit(2, "line"),
legend.spacing.x = unit(0.1, 'cm'),
legend.position = "bottom",
)
}
# Make default base map plot
xmin2 <- 303000; xmax2 <- 916000; xrange <- xmax2 - xmin2
ymin2 <- 5980000; ymax2 <- 6450000; yrange <- ymax2 - ymin2
plot_map_raster <-
ggplot(swe_coast_proj) +
xlim(xmin2, xmax2) +
ylim(ymin2, ymax2) +
labs(x = "Longitude", y = "Latitude") +
geom_sf(size = 0.3) +
theme_plot() +
guides(colour = guide_colorbar(title.position = "top", title.hjust = 0.5),
fill = guide_colorbar(title.position = "top", title.hjust = 0.5))
# Read the split files and join them
d1 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod-condition/master/data/for_analysis/mdat_cond_(1_2).csv")
d2 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod-condition/master/data/for_analysis/mdat_cond_(2_2).csv")
d <- bind_rows(d1, d2)
unique(is.na(d$density_cod))
#> [1] FALSE
unique(is.na(d$density_cod_rec))
#> [1] FALSE TRUE
# Calculate standardized variables
d <- d %>%
mutate(ln_length_cm = log(length_cm),
ln_weight_g = log(weight_g),
year_ct = year - mean(year),
biomass_her_sc = biomass_her,
biomass_her_sd_sc = biomass_her_sd,
biomass_spr_sc = biomass_spr,
biomass_spr_sd_sc = biomass_spr_sd,
density_cod_sc = density_cod,
density_cod_rec_sc = density_cod_rec,
density_fle_sc = density_fle,
density_fle_rec_sc = density_fle_rec,
density_saduria_sc = density_saduria,
density_saduria_rec_sc = density_saduria_rec,
depth_sc = depth,
depth_rec_sc = depth_rec,
oxy_sc = oxy,
oxy_rec_sc = oxy_rec,
temp_sc = temp,
temp_rec_sc = temp_rec,
year_f = as.factor(year)) %>%
mutate_at(c("biomass_her_sc", "biomass_her_sd_sc", "biomass_spr_sc", "biomass_spr_sd_sc",
"density_cod_sc", "density_cod_rec_sc",
"density_fle_sc", "density_fle_rec_sc",
"density_saduria_sc", "density_saduria_rec_sc",
"depth_sc", "depth_rec_sc",
"oxy_sc", "oxy_rec_sc",
"temp_sc", "temp_rec_sc"
),
~(scale(.) %>% as.vector)) %>%
mutate(year = as.integer(year))
unique(is.na(d))
#> weight_g length_cm year lat lon ices_rect sub_div depth oxy temp
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [3,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [4,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [5,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [6,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [7,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> X Y density_cod density_fle density_cod_data density_fle_data
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [3,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [4,] FALSE FALSE FALSE FALSE TRUE TRUE
#> [5,] FALSE FALSE FALSE FALSE TRUE TRUE
#> [6,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [7,] FALSE FALSE FALSE FALSE TRUE TRUE
#> density_saduria density_saduria_rec density_cod_rec density_fle_rec
#> [1,] FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE
#> [3,] FALSE FALSE FALSE FALSE
#> [4,] FALSE FALSE FALSE FALSE
#> [5,] FALSE TRUE TRUE TRUE
#> [6,] FALSE TRUE TRUE TRUE
#> [7,] FALSE FALSE FALSE FALSE
#> depth_rec oxy_rec temp_rec biomass_spr biomass_her biomass_spr_sd
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE TRUE TRUE FALSE
#> [3,] FALSE FALSE FALSE TRUE TRUE TRUE
#> [4,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [5,] TRUE TRUE TRUE TRUE TRUE FALSE
#> [6,] TRUE TRUE TRUE TRUE TRUE FALSE
#> [7,] FALSE FALSE FALSE TRUE TRUE FALSE
#> biomass_her_sd ln_length_cm ln_weight_g year_ct biomass_her_sc
#> [1,] FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE TRUE
#> [3,] TRUE FALSE FALSE FALSE TRUE
#> [4,] FALSE FALSE FALSE FALSE FALSE
#> [5,] FALSE FALSE FALSE FALSE TRUE
#> [6,] FALSE FALSE FALSE FALSE TRUE
#> [7,] FALSE FALSE FALSE FALSE TRUE
#> biomass_her_sd_sc biomass_spr_sc biomass_spr_sd_sc density_cod_sc
#> [1,] FALSE FALSE FALSE FALSE
#> [2,] FALSE TRUE FALSE FALSE
#> [3,] TRUE TRUE TRUE FALSE
#> [4,] FALSE FALSE FALSE FALSE
#> [5,] FALSE TRUE FALSE FALSE
#> [6,] FALSE TRUE FALSE FALSE
#> [7,] FALSE TRUE FALSE FALSE
#> density_cod_rec_sc density_fle_sc density_fle_rec_sc density_saduria_sc
#> [1,] FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE
#> [3,] FALSE FALSE FALSE FALSE
#> [4,] FALSE FALSE FALSE FALSE
#> [5,] TRUE FALSE TRUE FALSE
#> [6,] TRUE FALSE TRUE FALSE
#> [7,] FALSE FALSE FALSE FALSE
#> density_saduria_rec_sc depth_sc depth_rec_sc oxy_sc oxy_rec_sc temp_sc
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [3,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [4,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [5,] TRUE FALSE TRUE FALSE TRUE FALSE
#> [6,] TRUE FALSE TRUE FALSE TRUE FALSE
#> [7,] FALSE FALSE FALSE FALSE FALSE FALSE
#> temp_rec_sc year_f
#> [1,] FALSE FALSE
#> [2,] FALSE FALSE
#> [3,] FALSE FALSE
#> [4,] FALSE FALSE
#> [5,] TRUE FALSE
#> [6,] TRUE FALSE
#> [7,] FALSE FALSE
unique(is.na(d$density_cod_rec))
#> [1] FALSE TRUE
unique(is.na(d$density_cod))
#> [1] FALSE
# Drop NA covariates for the correlation plot
d <- d %>% drop_na(biomass_her_sc, biomass_her_sd_sc, biomass_spr_sc, biomass_spr_sd_sc,
density_cod_sc, density_cod_rec_sc, density_fle_sc, density_fle_rec_sc,
density_saduria_sc, density_saduria_rec_sc, depth_sc, depth_rec_sc,
oxy_sc, oxy_rec_sc, temp_sc, temp_rec_sc)
# Sample size
nrow(d)
#> [1] 94295
# stdev of oxygen
sd(d$oxy)
#> [1] 1.805214
# Plot correlation between variables
d_cor <- d %>% dplyr::select("biomass_her_sc", "biomass_her_sd_sc",
"biomass_spr_sc", "biomass_spr_sd_sc",
"density_cod_sc", "density_cod_rec_sc",
"density_fle_sc", "density_fle_rec_sc",
"density_saduria_sc", "density_saduria_rec_sc",
"depth_sc", "depth_rec_sc",
"oxy_sc", "oxy_rec_sc",
"temp_sc", "temp_rec_sc")
corr <- round(cor(d_cor), 1)
ggcorrplot(corr, hc.order = TRUE, type = "lower", lab = TRUE, lab_size = 2.5) +
theme_classic() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.3))
ggsave("figures/supp/condition/correlation_variables.png", width = 6.5, height = 6.5, dpi = 600)
# Plot data
d_plot <- d %>%
group_by(year, X, Y) %>%
summarise(n = n())
plot_map_raster +
geom_point(data = d_plot, aes(X*1000, Y*1000, size = n, color = n), alpha = 0.5) +
scale_color_viridis() +
scale_size(range = c(.1, 3), name = "N per haul") +
facet_wrap(~year, ncol = 5) +
guides(size = "none") +
theme_facet_map()
ggsave("figures/supp/condition/spatial_cond_data.png", width = 6.5, height = 8.5, dpi = 600)
pred_grid1 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(1_2).csv")
pred_grid2 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(2_2).csv")
pred_grid <- bind_rows(pred_grid1, pred_grid2)
unique(is.na(pred_grid$density_cod))
#> [1] FALSE
unique(is.na(pred_grid$density_cod_rec))
#> [1] FALSE
pred_grid <- pred_grid %>% mutate(year = as.integer(year),
year_f = as.factor(year))
# Scale the variables with respect to data! First calculate means in data
data_means <- d %>%
dplyr::select(biomass_her, biomass_her_sd, biomass_spr, biomass_spr_sd,
density_cod, density_cod_rec,
density_fle, density_fle_rec,
density_saduria, density_saduria_rec,
depth, depth_rec,
oxy, oxy_rec,
temp, temp_rec) %>%
mutate_all(~mean(.)) %>%
distinct(.keep_all = TRUE)
data_stdev <- d %>%
dplyr::select(biomass_her, biomass_her_sd, biomass_spr, biomass_spr_sd,
density_cod, density_cod_rec,
density_fle, density_fle_rec,
density_saduria, density_saduria_rec,
depth, depth_rec,
oxy, oxy_rec,
temp, temp_rec
) %>%
mutate_all(~sd(.)) %>%
distinct(.keep_all = TRUE)
# Before actually scaling, replace the ices-rectangle pelagic values that are NA with the mean across rectangles in the sub division.
# Replace the sub-division pelagic values that are NA with the mean in the year.
# Note that the sub-division values are not the sum of the rectangles (due to missing rectangles), so
# I need to calculate a sub-division mean across rectangles within each sub-division
pred_grid <- pred_grid %>%
group_by(year) %>%
mutate(median_biomass_sprat_across_sub_div = median(biomass_spr_sd, na.rm = TRUE),
median_biomass_herring_across_sub_div = median(biomass_her_sd, na.rm = TRUE)) %>%
ungroup() %>% # Replace sub_divsion NA's with the median across sub_divisions in that year
mutate(biomass_spr_sd = ifelse(is.na(biomass_spr_sd) == TRUE, median_biomass_sprat_across_sub_div, biomass_spr_sd),
biomass_her_sd = ifelse(is.na(biomass_her_sd) == TRUE, median_biomass_herring_across_sub_div, biomass_her_sd)) %>%
group_by(year, sub_div) %>%
mutate(median_biomass_sprat_across_rect_in_sub_div = median(biomass_spr, na.rm = TRUE),
median_biomass_herring_across_rect_in_sub_div = median(biomass_her, na.rm = TRUE)) %>%
ungroup() %>% # Replace rectangle NA's with the median across rectangles in that year and sub-division
mutate(biomass_spr = ifelse(is.na(biomass_spr) == TRUE, median_biomass_sprat_across_rect_in_sub_div, biomass_spr),
biomass_her = ifelse(is.na(biomass_her) == TRUE, median_biomass_herring_across_rect_in_sub_div, biomass_her)) %>%
# Since I still have some NAs (some sub-divisions do not have a single rectangle in some years), I will will those rectangles
# with the sub-division value divided by the number of rectangles in that sub division
group_by(year, sub_div) %>%
mutate(n_rect = length(unique(ices_rect))) %>%
ungroup() %>%
mutate(biomass_spr = ifelse(is.na(biomass_spr) == TRUE, biomass_spr_sd/n_rect, biomass_spr),
biomass_her = ifelse(is.na(biomass_her) == TRUE, biomass_her_sd/n_rect, biomass_her))
pred_grid <- pred_grid %>%
mutate(ln_length_cm = log(1),
year_ct = year - mean(year),
biomass_her_sc = (biomass_her - data_means$biomass_her) / data_stdev$biomass_her,
biomass_her_sd_sc = (biomass_her_sd - data_means$biomass_her_sd) / data_stdev$biomass_her_sd,
biomass_spr_sc = (biomass_spr - data_means$biomass_spr) / data_stdev$biomass_spr,
biomass_spr_sd_sc = (biomass_spr_sd - data_means$biomass_spr_sd) / data_stdev$biomass_spr_sd,
density_cod_sc = (density_cod - data_means$density_cod) / data_stdev$density_cod,
density_cod_rec_sc = (density_cod_rec - data_means$density_cod_rec) / data_stdev$density_cod_rec,
density_fle_sc = (density_fle - data_means$density_fle) / data_stdev$density_fle,
density_fle_rec_sc = (density_fle_rec - data_means$density_fle_rec) / data_stdev$density_fle_rec,
density_saduria_sc = (density_saduria - data_means$density_saduria) / data_stdev$density_saduria,
density_saduria_rec_sc = (density_saduria_rec - data_means$density_saduria_rec) / data_stdev$density_saduria_rec,
depth_sc = (depth - data_means$depth) / data_stdev$depth,
depth_rec_sc = (depth_rec - data_means$depth_rec) / data_stdev$depth_rec,
oxy_sc = (oxy - data_means$oxy) / data_stdev$oxy,
oxy_rec_sc = (oxy_rec - data_means$oxy_rec) / data_stdev$oxy_rec,
temp_sc = (temp - data_means$temp) / data_stdev$temp,
temp_rec_sc = (temp_rec - data_means$temp_rec) / data_stdev$temp_rec,
)
# Plot on map:
ggplot(pred_grid2, aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)
pred_grid %>%
drop_na() %>%
ggplot(., aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)
pred_grid %>%
drop_na(biomass_spr_sd) %>%
ggplot(., aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)
unique(is.na(pred_grid))
#> X Y depth year id oxy_q3 oxy temp_q3 temp lon lat
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> sub_div ices_rect density_saduria biomass_spr biomass_her biomass_spr_sd
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE FALSE
#> biomass_her_sd density_cod density_fle depth_rec temp_rec oxy_rec
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE FALSE
#> density_cod_rec density_fle_rec density_saduria_rec year_f
#> [1,] FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE
#> median_biomass_sprat_across_sub_div median_biomass_herring_across_sub_div
#> [1,] FALSE FALSE
#> [2,] FALSE FALSE
#> median_biomass_sprat_across_rect_in_sub_div
#> [1,] FALSE
#> [2,] TRUE
#> median_biomass_herring_across_rect_in_sub_div n_rect ln_length_cm year_ct
#> [1,] FALSE FALSE FALSE FALSE
#> [2,] TRUE FALSE FALSE FALSE
#> biomass_her_sc biomass_her_sd_sc biomass_spr_sc biomass_spr_sd_sc
#> [1,] FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE
#> density_cod_sc density_cod_rec_sc density_fle_sc density_fle_rec_sc
#> [1,] FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE
#> density_saduria_sc density_saduria_rec_sc depth_sc depth_rec_sc oxy_sc
#> [1,] FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE
#> oxy_rec_sc temp_sc temp_rec_sc
#> [1,] FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE
pred_grid %>%
drop_na(biomass_spr_sc, biomass_spr_sd_sc, biomass_her_sc, biomass_her_sd_sc) %>%
ggplot(., aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)
spde <- make_mesh(d, xy_cols = c("X", "Y"),
n_knots = 100,
type = "kmeans", seed = 42)
# Plot and save spde
png(file = "figures/supp/condition/spde.png", units = "in", width = 6.5, height = 6.5, res = 300)
plot(spde)
dev.off()
ptm <- proc.time()
mnull <- sdmTMB(
formula = ln_weight_g ~ ln_length_cm,
data = d,
time = NULL,
mesh = spde,
family = student(link = "identity", df = 5),
spatiotemporal = "off",
spatial = "off",
silent = TRUE,
reml = FALSE
)
proc.time() - ptm
#> user system elapsed
#> 1.611 0.116 1.763
sanity(mnull)
#> ✖ Non-linear minimizer did not converge: do not trust this model!
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
# mlm <- lm(
# formula = ln_weight_g ~ ln_length_cm,
# data = d)
# bm <- brms::brm(formula = ln_weight_g ~ ln_length_cm, data = d, family = student())
# summary(bm)
# plot(bm)
# brms::pp_check(bm)
# Right, so even though sdmTMB warns, the actual estimates are very similar to lm and brms, student-t or gaussian. So will just play along with it
# Extract average a and b across the total dataset
a <- as.numeric(tidy(mnull)[1, "estimate"])
b <- as.numeric(tidy(mnull)[2, "estimate"])
a
#> [1] -4.638609
b
#> [1] 2.984044
d$weight_g_avg_pred <- exp(a)*d$length_cm^b
plot(d$weight_g ~ d$length_cm)
plot(d$weight_g_avg_pred ~ d$length_cm)
plot(d$weight_g_avg_pred ~ d$weight_g); abline(a = 0, b = 1, col = "red")
d$le_cren <- d$weight_g / d$weight_g_avg_pred
d$log_le_cren <- log(d$le_cren)
hist(d$le_cren)
ptm <- proc.time()
mfull <- sdmTMB(
formula = log_le_cren ~ -1 + year_f + biomass_her_sc + biomass_her_sd_sc +
biomass_spr_sc + biomass_spr_sd_sc +
density_cod_sc + density_cod_rec_sc +
density_fle_sc + density_fle_rec_sc +
density_saduria_sc + density_saduria_rec_sc +
depth_sc + depth_rec_sc +
oxy_sc + oxy_rec_sc +
temp_sc + temp_rec_sc,
data = d,
time = "year",
mesh = spde,
family = student(link = "identity", df = 5),
spatiotemporal = "iid",
spatial = "on",
silent = TRUE,
reml = FALSE,
control = sdmTMBcontrol(newton_loops = 1)#,
# priors = sdmTMBpriors(b = normal(c(rep(NA, 27), rep(0, 16)),
# c(rep(NA, 27), rep(1, 16))))
)
system("say Just finished!")
proc.time() - ptm
#> user system elapsed
#> 883.707 43.842 944.490
sanity(mfull)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mfull, conf.int = TRUE) %>% arrange(desc(abs(estimate)))
tidy(mfull, effects = "ran_par", conf.int = TRUE)
# Arrange coefficients
tidy(mfull, conf.int = TRUE) %>% filter(!grepl("year", term)) %>% arrange(desc(estimate))
#> filter: removed 27 rows (63%), 16 rows remaining
tidy(mfull, conf.int = TRUE) %>% filter(!grepl("year", term)) %>% arrange(estimate)
#> filter: removed 27 rows (63%), 16 rows remaining
# Average magnitude of fixed effects?
tidy(mfull, conf.int = TRUE) %>%
filter(!grepl("year", term)) %>%
mutate(abs_coef = abs(estimate)) %>%
summarise(mean_abs_coef = mean(abs_coef))
#> filter: removed 27 rows (63%), 16 rows remaining
#> mutate: new variable 'abs_coef' (double) with 16 unique values and 0% NA
#> summarise: now one row and one column, ungrouped
tidy(mfull, conf.int = TRUE, effects = "ran_par") %>%
filter(term %in% c("sigma_O", "sigma_E")) %>%
summarise(mean_estimate = mean(estimate))
#> filter: removed 2 rows (50%), 2 rows remaining
#> summarise: now one row and one column, ungrouped
# This model is for the sensitivity analysis, comparing estimates between different model subsets etc. These tend to fit better with only a spatiotemporal term (spatial stdev tends to blow up), so I want them to be as comparable as possible
ptm <- proc.time()
mfull_spatiotemporal <- sdmTMB(
formula = log_le_cren ~ -1 + year_f + biomass_her_sc + biomass_her_sd_sc +
biomass_spr_sc + biomass_spr_sd_sc +
density_cod_sc + density_cod_rec_sc +
density_fle_sc + density_fle_rec_sc +
density_saduria_sc + density_saduria_rec_sc +
depth_sc + depth_rec_sc +
oxy_sc + oxy_rec_sc +
temp_sc + temp_rec_sc,
data = d,
time = "year",
mesh = spde,
family = student(link = "identity", df = 5),
spatiotemporal = "iid",
spatial = "off",
silent = TRUE,
reml = FALSE,
control = sdmTMBcontrol(newton_loops = 1)
)
system("say Just finished!")
proc.time() - ptm
#> user system elapsed
#> 672.519 20.072 704.688
sanity(mfull_spatiotemporal)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mfull_spatiotemporal, conf.int = TRUE) %>% arrange(desc(abs(estimate)))
tidy(mfull_spatiotemporal, effects = "ran_par", conf.int = TRUE)
# Extract random and fixed coefficients from the full model
mfull_est_st <- bind_rows(tidy(mfull_spatiotemporal, effects = "ran_par", conf.int = TRUE) %>%
filter(term == "sigma_E"),
tidy(mfull_spatiotemporal, effects = "fixed", conf.int = TRUE) %>%
filter(!grepl('year', term))) %>%
mutate(term = factor(term))
#> filter: removed 2 rows (67%), one row remaining
#> filter: removed 27 rows (63%), 16 rows remaining
#> mutate: converted 'term' from character to factor (0 new NA)
# Extract coefficients and save as csv
mfull_est_st <- mfull_est_st %>%
mutate(group = "Herring",
group = ifelse(term %in% c("biomass_spr_sc", "biomass_spr_sd_sc"), "Sprat", group),
group = ifelse(term %in% c("density_cod_rec_sc", "density_cod_sc"), "Cod", group),
group = ifelse(term %in% c("density_fle_rec_sc", "density_fle_sc"), "Flounder", group),
group = ifelse(term %in% c("density_saduria_rec_sc", "density_saduria_sc"), "Saduria", group),
group = ifelse(term %in% c("depth_rec_sc", "depth_sc"), "Depth", group),
group = ifelse(term %in% c("oxy_rec_sc", "oxy_sc"), "Oxygen", group),
group = ifelse(term %in% c("temp_rec_sc", "temp_sc"), "Temp", group),
group = ifelse(term == "sigma_E", "Random", group),
) %>%
mutate(scale = "Subdivision",
scale = ifelse(term %in% c("biomass_her_sc", "biomass_spr_sc", "density_cod_rec_sc", "density_fle_rec_sc",
"density_saduria_rec_sc", "depth_rec_sc", "oxy_rec_sc", "temp_rec_sc"), "Rectangle", scale),
scale = ifelse(term %in% c("density_cod_sc", "density_fle_sc", "density_saduria_sc",
"depth_sc", "oxy_sc", "temp_sc"), "Haul", scale),
scale = ifelse(term == "sigma_E", "Spatial/spatiotemporal s.d.", scale),
) %>%
mutate(term2 = ifelse(term == "sigma_E", 2, 1))
#> mutate: new variable 'group' (character) with 9 unique values and 0% NA
#> mutate: new variable 'scale' (character) with 4 unique values and 0% NA
#> mutate: new variable 'term2' (double) with 2 unique values and 0% NA
# Now save this for comparison in sensitivity analysis script
write.csv(mfull_est_st, "output/mfull_est_st.csv")
# Save model (so that I can add density predictions to the pred grid for condition)
saveRDS(mfull, "output/mfull.rds")
mcmc_res <- residuals(mfull, type = "mle-mcmc", mcmc_iter = 5000, mcmc_warmup = 2000)
png(file = "figures/supp/condition/qq.png", units = "in", width = 6.5, height = 6.5, res = 300)
qqnorm(mcmc_res);qqline(mcmc_res)
dev.off()
d$residuals <- mcmc_res[, 1]
r2.sdmTMBb_fe <- tidy(mfull)
b_fe <- stats::setNames(b_fe$estimate, b_fe$term)
X <- mfull$tmb_data$X_ij
X <- matrix(unlist(X), ncol = length(b_fe))
VarF <- var(as.vector(b_fe %*% t(X))) # variance from fixed-effects
# What if I skip the year effect?
b_fe2 <- tidy(mfull) %>% filter(!grepl('year', term))
#> filter: removed 27 rows (63%), 16 rows remaining
b_fe2 <- stats::setNames(b_fe2$estimate, b_fe2$term)
X2 <- mfull$tmb_data$X_ij
first_index <- match('year_f2019', names(as.data.frame(mfull$tmb_data$X_ij[[1]]))) + 1
last_index <- ncol(mfull$tmb_data$X_ij[[1]])
X2 <- matrix(unlist(X2), ncol = length(b_fe))[, first_index:last_index] # skip year columns
VarF2 <- var(as.vector(b_fe2 %*% t(X2))) # variance from fixed-effects
b_ran <- tidy(mfull, "ran_par")
sigma_O <- b_ran$estimate[b_ran$term == "sigma_O"] # spatial standard deviation
sigma_E <- b_ran$estimate[b_ran$term == "sigma_E"] # spatiotemporal standard deviation
VarO <- sigma_O^2 # spatial variance
VarE <- sigma_E^2 # spatiotemporal variance
# Calculate obs - pred
d_r2 <- d
d_r2$pred <- predict(mfull)
d_r2$resid_manual <- d_r2$le_cren - d_r2$pred$est
#hist(d_r2$resid_manual)
VarR <- var(as.vector(d_r2$resid_manual)) # "residual" variance
# Marginal R2s
denominator <- VarF + VarO + VarE + VarR
denominator_no_yr <- VarF2 + VarO + VarE + VarR
marg <- VarF/denominator
marg_no_yr <- VarF2/denominator_no_yr
out <- data.frame(
marginal = marg,
partial_spatial = VarO/denominator,
partial_spatiotemporal = VarE/denominator,
conditional = (VarF + VarO + VarE)/denominator,
all_random = (VarO + VarE)/denominator,
marginal_random_ratio = marg / ((VarO + VarE)/denominator),
random_marginal_ratio = ((VarO + VarE)/denominator) / marg
)
out
out_no_yr <- data.frame(
marginal = marg_no_yr,
partial_spatial = VarO/denominator_no_yr,
partial_spatiotemporal = VarE/denominator_no_yr,
conditional = (VarF2 + VarO + VarE)/denominator_no_yr,
all_random = (VarO + VarE)/denominator_no_yr,
marginal_random_ratio = marg_no_yr / ((VarO + VarE)/denominator_no_yr),
random_marginal_ratio = ((VarO + VarE)/denominator_no_yr) / marg_no_yr
)
out_no_yr
# Residuals on map
plot_map_raster +
geom_point(data = d, aes(x = X * 1000, y = Y * 1000, color = residuals), size = 0.5) +
scale_colour_gradient2() +
facet_wrap(~ year, ncol = 5) +
labs(color = "residuals") +
theme_facet_map() +
geom_sf(size = 0.1)
ggsave("figures/supp/condition/spatial_resid.png", width = 6.5, height = 8.5, dpi = 600)
# Residuals vs length and year
ggplot(d, aes(x = length_cm, y = residuals, color = residuals)) +
geom_point(size = 0.1, alpha = 1) +
geom_hline(yintercept = 0, size = 0.5, color = "black", linetype = 2, alpha = 0.5) +
labs(x = "length [cm]", y = "Residuals", color = "") +
stat_smooth(color = "black", size = 0.5) +
facet_wrap(~year, ncol = 5) +
scale_color_gradient2() +
theme_plot()
ggsave("figures/supp/condition/residuals_vs_length_and_year.png", width = 6.5, height = 8.5, dpi = 600)
ggplot(d, aes(year, length_cm, z = residuals)) +
stat_summary_2d(bins = 40) +
scale_fill_gradient2()
# ggsave("figures/supp/condition/residuals_vs_length_and_year2.png", width = 6.5, height = 6.5, dpi = 600)
get_index_simspred_avg_sim <- predict(mfull, newdata = pred_grid %>% filter(depth > 10 & depth < 110), nsim = 500)
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
# Weighted version. We could do it directly in the get_index_sims or manually in the prediction grid (but then we don't get an CIs). This is so that we can divide by the sum of weights (needs to be done by year)
weight_sum <- pred_grid %>%
filter(depth > 10 & depth < 110) %>%
group_by(year) %>%
summarise(density_cod = sum(density_cod))
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 2 columns, ungrouped
index_avg_sim <- get_index_sims(pred_avg_sim,
area = filter(pred_grid, depth > 10 & depth < 110)$density_cod) %>%
left_join(weight_sum) %>%
mutate(
est = est / density_cod,
lwr = lwr / density_cod,
upr = upr / density_cod
)
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Joining, by = "year"
#> left_join: added one column (density_cod)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 27
#> > ====
#> > rows total 27
#> mutate: changed 27 values (100%) of 'est' (0 new NA)
#> changed 27 values (100%) of 'lwr' (0 new NA)
#> changed 27 values (100%) of 'upr' (0 new NA)
index_avg_sims <- get_index_sims(pred_avg_sim,
area = filter(pred_grid, depth > 10 & depth < 110)$density_cod,
return_sims = TRUE) %>%
left_join(weight_sum) %>%
mutate(est = .value / density_cod)
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Joining, by = "year"
#> left_join: added one column (density_cod)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 13,500
#> > ========
#> > rows total 13,500
#> mutate: new variable 'est' (double) with 13,500 unique values and 0% NA
ggplot(index_avg_sim, aes(year, est, ymin = lwr, ymax = upr)) +
geom_ribbon(alpha = 0.2, color = NA) +
geom_line(size = 1) +
geom_line(data = filter(index_avg_sims, .iteration < 26),
aes(year, est, group = .iteration), inherit.aes = FALSE, alpha = 0.3)
#> filter: removed 12,825 rows (95%), 675 rows remaining
# Make the un-weighted index
ncells <- filter(pred_grid, year == max(pred_grid$year) & depth > 10 & depth < 110) %>% nrow()
#> filter: removed 250,890 rows (97%), 7,689 rows remaining
index_avg_sim_uw <- get_index_sims(pred_avg_sim, area = rep(1/ncells, nrow(pred_avg_sim)))
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
index_avg_sims_uw <- get_index_sims(pred_avg_sim, area = rep(1/ncells, nrow(pred_avg_sim)),
return_sims = TRUE) # This is for calculation on samples!
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
pred_grid_fe <- pred_grid %>%
mutate(year_f = as.factor(1993))
#> mutate: changed 249,002 values (96%) of 'year_f' (0 new NA)
pred_avg_sim_fe <- predict(mfull, newdata = pred_grid_fe %>% filter(depth > 10 & depth < 110), nsim = 500, re_form = NA)
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
weight_sum <- pred_grid %>%
filter(depth > 10 & depth < 110) %>%
group_by(year) %>%
summarise(density_cod = sum(density_cod))
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 2 columns, ungrouped
index_avg_sim_fe <- get_index_sims(pred_avg_sim_fe,
area = filter(pred_grid, depth > 10 & depth < 110)$density_cod) %>%
left_join(weight_sum) %>%
mutate(
est = est / density_cod,
lwr = lwr / density_cod,
upr = upr / density_cod
)
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Joining, by = "year"
#> left_join: added one column (density_cod)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 27
#> > ====
#> > rows total 27
#> mutate: changed 27 values (100%) of 'est' (0 new NA)
#> changed 27 values (100%) of 'lwr' (0 new NA)
#> changed 27 values (100%) of 'upr' (0 new NA)
index_avg_sims_fe <- get_index_sims(pred_avg_sim_fe,
area = filter(pred_grid, depth > 10 & depth < 110)$density_cod,
return_sims = TRUE) %>%
left_join(weight_sum) %>%
mutate(est = .value / density_cod)
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Joining, by = "year"
#> left_join: added one column (density_cod)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 13,500
#> > ========
#> > rows total 13,500
#> mutate: new variable 'est' (double) with 13,500 unique values and 0% NA
# With random effects but no fixed year effect
# pred_grid_re <- pred_grid %>%
# mutate(year_f = as.factor(1993))
#
# pred_avg_sim_re <- predict(mfull, newdata = pred_grid_re %>% filter(depth > 10 & depth < 110), nsim = 500)
#
# weight_sum <- pred_grid %>%
# filter(depth > 10 & depth < 110) %>%
# group_by(year) %>%
# summarise(density_cod = sum(density_cod))
#
# index_avg_sim_re <- get_index_sims(pred_avg_sim_re,
# area = filter(pred_grid, depth > 10 & depth < 110)$density_cod) %>%
# left_join(weight_sum) %>%
# mutate(
# est = est / density_cod,
# lwr = lwr / density_cod,
# upr = upr / density_cod
# )
#
# index_avg_sims_re <- get_index_sims(pred_avg_sim_re,
# area = filter(pred_grid, depth > 10 & depth < 110)$density_cod,
# return_sims = TRUE) %>%
# left_join(weight_sum) %>%
# mutate(est = .value / density_cod)
# Plot all together
index_comp_sim <- bind_rows(#index_avg_sim_re %>% mutate(Prediction = "Random & Fixed effects (no yr)"),
index_avg_sim_fe %>% mutate(Prediction = "Fixed effects (no yr)"),
index_avg_sim %>% mutate(Prediction = "Full"))
#> mutate: new variable 'Prediction' (character) with one unique value and 0% NA
#> mutate: new variable 'Prediction' (character) with one unique value and 0% NA
index_comp_sims <- bind_rows(#index_avg_sims_re %>% mutate(Prediction = "Random & Fixed effects (no yr)"),
index_avg_sims_fe %>% mutate(Prediction = "Fixed effects (no yr)"),
index_avg_sims %>% mutate(Prediction = "Full")) %>%
filter(.iteration < 26)
#> mutate: new variable 'Prediction' (character) with one unique value and 0% NA
#> mutate: new variable 'Prediction' (character) with one unique value and 0% NA
#> filter: removed 25,650 rows (95%), 1,350 rows remaining
ggplot(index_comp_sim, aes(year, est, ymin = lwr, ymax = upr, color = Prediction, fill = Prediction)) +
geom_ribbon(alpha = 0.2, color = NA) +
geom_line(size = 1) +
geom_line(data = filter(index_comp_sims, .iteration < 26),
aes(year, est, group = interaction(.iteration, Prediction), color = Prediction), inherit.aes = FALSE, alpha = 0.3) +
NULL
#> filter: no rows removed
# Reshape data again to calculate row-wise differences
year_diff <- index_avg_sims %>%
dplyr::select(-.value, -density_cod) %>%
filter(year %in% c(1993, 2019)) %>%
pivot_wider(names_from = year, values_from = est) %>%
mutate("2019-1993" = `2019` - `1993`) %>%
pivot_longer(2:4) %>%
ungroup()
#> filter: removed 12,500 rows (93%), 1,000 rows remaining
#> pivot_wider: reorganized (year, est) into (1993, 2019) [was 1000x3, now 500x3]
#> mutate: new variable '2019-1993' (double) with 500 unique values and 0% NA
#> pivot_longer: reorganized (1993, 2019, 2019-1993) into (name, value) [was 500x4, now 1500x3]
#> ungroup: no grouping variables
# Plot these
p1 <- ggplot(year_diff) +
geom_density(data = filter(year_diff, !name == "2019-1993"),
aes(value, fill = name), alpha = 0.6) +
scale_fill_brewer(palette = "Dark2") +
coord_cartesian(expand = 0)
#> filter: removed 500 rows (33%), 1,000 rows remaining
p2 <- ggplot(year_diff) +
geom_density(data = filter(year_diff, name == "2019-1993"),
aes(value), alpha = 0.6, fill = "grey") +
geom_vline(xintercept = 0, linetype = 2)
#> filter: removed 1,000 rows (67%), 500 rows remaining
p1 / p2
# Calculate quantiles for each level (to go in the results section)
year_diff %>%
group_by(name) %>%
summarise(median = quantile(value, probs = 0.5),
upr97.5 = quantile(value, probs = 0.975),
lwr2.5 = quantile(value, probs = 0.025))
#> group_by: one grouping variable (name)
#> summarise: now 3 rows and 4 columns, ungrouped
# Percent change in condition factor
year_diff %>%
group_by(name) %>%
pivot_wider(names_from = name, values_from = value) %>%
mutate(perc_change = ((`2019`-`1993`)/`1993`)*100) %>%
ungroup() %>%
summarise(lwr2.5 = quantile(perc_change, probs = 0.025),
median = quantile(perc_change, probs = 0.5),
upr97.5 = quantile(perc_change, probs = 0.975))
#> group_by: one grouping variable (name)
#> pivot_wider: reorganized (name, value) into (1993, 2019, 2019-1993) [was 1500x3, now 500x4]
#> mutate: new variable 'perc_change' (double) with 500 unique values and 0% NA
#> ungroup: no grouping variables
#> summarise: now one row and 3 columns, ungrouped
Evaluate sensitivity to using a pred grid with a subset excluding low density areas
# Now do the same but excluding the areas not sampled in the pred grid
# ncells_sub <- filter(pred_grid, year == max(pred_grid$year) & depth > 40 & depth < 120) %>% nrow()
# pred_grid_sub <- filter(pred_grid, depth > 40 & depth < 120)
#
# pred_avg_sim_sub <- predict(mfull, newdata = pred_grid_sub, nsim = 100)
#
# index_avg_sim_sub <- get_index_sims(pred_avg_sim_sub, area = rep(1/ncells_sub, nrow(pred_avg_sim_sub)))
#
# # Combine and plot & compare
# index_avg_sim_comp <- bind_rows(mutate(index_avg_sim, area = "full"),
# mutate(index_avg_sim_sub, area = "subset"))
#
# ggplot(index_avg_sim_comp, aes(year, est, ymin = lwr, ymax = upr, color = area, fill = area)) +
# geom_line() +
# geom_ribbon(alpha = 0.4, color = NA)
div_index_list <- list()
for(i in unique(pred_grid$sub_div)){
pred_grid_sub <- pred_grid %>% filter(sub_div == i & depth > 10 & depth < 110)
pred_avg_sim_sub <- predict(mfull, newdata = pred_grid_sub, nsim = 500)
weight_sum_sub <- pred_grid_sub %>%
group_by(year) %>%
summarise(density_cod = sum(density_cod))
index_avg_sim <- get_index_sims(pred_avg_sim_sub,
area = pred_grid_sub$density_cod) %>%
left_join(weight_sum_sub) %>%
mutate(est = est / density_cod,
lwr = lwr / density_cod,
upr = upr / density_cod) %>%
mutate(sub_div = i)
div_index_list[[i]] <- index_avg_sim
}
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
div_index <- bind_rows(div_index_list)
# Spatiotemporally standardized indicies
div_index2 <- bind_rows(div_index %>% dplyr::select(year, est, lwr, upr, sub_div) %>% mutate(sub_div = as.factor(sub_div)),
mutate(index_avg_sim, sub_div = as.factor("Weighted index")))
# Add in raw mean of data:
div_index2 <- bind_rows(div_index2,
d %>%
group_by(year) %>%
summarise(est = exp(mean(log_le_cren))) %>%
mutate(sub_div = "Empirical mean"))
# Add fixed effect predictions and the un-weighted index
div_index2 <- bind_rows(div_index2,
index_avg_sim_fe %>% mutate(sub_div = "Weighted index (fixed effects)"),
index_avg_sim_uw %>% mutate(sub_div = "Un-weighted index",
upr = NA,
lwr = NA))
# Add sims for total and fixed effect prediction
index_comp_sims <- bind_rows(index_avg_sims_fe %>% mutate(sub_div = "Weighted index (fixed effects)"),
index_avg_sims %>% mutate(sub_div = "Weighted index")) %>%
filter(.iteration < 26)
# Plot!
p1 <- div_index2 %>%
filter(sub_div %in% c("Weighted index", "Un-weighted index", "Empirical mean", "Weighted index (fixed effects)")) %>%
mutate(plot_facet = "Total and data") %>%
mutate(plot_group = ifelse(sub_div %in% c("Empirical mean", "Un-weighted index"), 1, 0)) %>%
mutate(line_group = ifelse(sub_div %in% c("Empirical mean", "Un-weighted index"), "1", "0")) %>%
ggplot(aes(year, est, color = reorder(sub_div, plot_group))) +
labs(y = "Le Cren's condition factor", x = "", color = "") +
geom_ribbon(aes(x = year, y = est, ymin = lwr, ymax = upr, fill = sub_div), color = NA, alpha = 0.25) +
geom_line(data = filter(index_comp_sims, .iteration < 26),
aes(year, est, group = interaction(.iteration, sub_div), color = sub_div), inherit.aes = FALSE, alpha = 0.3) +
geom_line(aes(linetype = line_group), size = 0.6) +
facet_wrap(~plot_facet, ncol = 1) +
guides(linetype = "none") +
scale_color_brewer(palette = "Set2") +
scale_fill_brewer(palette = "Set2") +
scale_x_continuous(breaks = c(1994, 2000, 2006, 2012, 2018)) +
theme_plot(base_size = 10) +
theme(plot.margin = unit(c(0.4, 0.4, 0, 0.4), "cm"),
legend.position = c(0.24, 0.11),
legend.direction = "horizontal",
legend.spacing.y = unit(0, 'cm'),
legend.key.size = unit(1, "cm"),
legend.key.width = unit(0.5, "cm"),
legend.title = element_text(size = 8),
legend.text = element_text(size = 7),
legend.background = element_rect(fill = NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()) +
guides(fill = "none",
color = guide_legend(ncol = 1, title.position = "top", title.hjust = 0.5))
p1
#> Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
#> -Inf
#> Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
#> -Inf
p2 <- div_index2 %>%
filter(!sub_div %in% c("Weighted index", "Un-weighted index", "Empirical mean", "Weighted index (fixed effects)")) %>%
mutate(plot_facet = "Subdivision") %>%
ggplot(aes(year, est, color = sub_div, fill = sub_div)) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
labs(y = "Le Cren's condition factor", x = "Year") +
geom_line() +
facet_wrap(~plot_facet, ncol = 1) +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2") +
scale_x_continuous(breaks = c(1994, 2000, 2006, 2012, 2018)) +
theme_plot(base_size = 10) +
theme(plot.margin = unit(c(0, 0.4, 0.4, 0.4), "cm"),
legend.position = c(0.92, 0.88),
legend.direction = "horizontal",
legend.spacing.y = unit(0, 'cm'),
legend.key.size = unit(0.8, "cm"),
legend.key.width = unit(0.3, "cm"),
legend.title = element_text(size = 8),
legend.text = element_text(size = 7),
legend.background = element_rect(fill = NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()) +
guides(color = guide_legend(ncol = 1, title.position = "top", title.hjust = 0.5),
fill = "none")
p2
p1 / p2 + plot_annotation(tag_levels = "A")
#> Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
#> -Inf
#> Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
#> -Inf
#ggsave("figures/cond_index.pdf", width = 11, height = 17, units = "cm")
# What is one standard deviation of oxygen?
sd(d$oxy)
#> [1] 1.805214
# Extract random and fixed coefficients from the full model
mfull_est <- bind_rows(tidy(mfull, effects = "ran_par", conf.int = TRUE) %>%
filter(term %in% c("sigma_O", "sigma_E")),
tidy(mfull, effects = "fixed", conf.int = TRUE) %>%
filter(!grepl('year', term))) %>%
mutate(term = factor(term))
mfull_est <- mfull_est %>%
mutate(group = "Herring",
group = ifelse(term %in% c("biomass_spr_sc", "biomass_spr_sd_sc"), "Sprat", group),
group = ifelse(term %in% c("density_cod_rec_sc", "density_cod_sc"), "Cod", group),
group = ifelse(term %in% c("density_fle_rec_sc", "density_fle_sc"), "Flounder", group),
group = ifelse(term %in% c("density_saduria_rec_sc", "density_saduria_sc"), "Saduria", group),
group = ifelse(term %in% c("depth_rec_sc", "depth_sc"), "Depth", group),
group = ifelse(term %in% c("oxy_rec_sc", "oxy_sc"), "Oxygen", group),
group = ifelse(term %in% c("temp_rec_sc", "temp_sc"), "Temp", group),
group = ifelse(term == "sigma_E", "Random", group),
group = ifelse(term == "sigma_O", "Random", group)
)
mfull_est <- mfull_est %>%
mutate(scale = "Subdivision",
scale = ifelse(term %in% c("biomass_her_sc", "biomass_spr_sc", "density_cod_rec_sc", "density_fle_rec_sc",
"density_saduria_rec_sc", "depth_rec_sc", "oxy_rec_sc", "temp_rec_sc"), "Rectangle", scale),
scale = ifelse(term %in% c("density_cod_sc", "density_fle_sc", "density_saduria_sc",
"depth_sc", "oxy_sc", "temp_sc"), "Haul", scale),
scale = ifelse(term == "sigma_E", "Spatial/spatiotemporal s.d.", scale),
scale = ifelse(term == "sigma_O", "Spatial/spatiotemporal s.d.", scale)
)
# Quick calculation:
mfull_est %>%
mutate(par_type = ifelse(term %in% c("sigma_O", "sigma_E"), "re", "fe")) %>%
group_by(par_type) %>%
summarise(mean_est = mean(abs(estimate))) %>%
pivot_wider(names_from = par_type, values_from = mean_est) %>%
mutate(ratio = re/fe)
# Plot effects
# This is the order changed labels should go in!
levels(mfull_est$term)
#> [1] "biomass_her_sc" "biomass_her_sd_sc" "biomass_spr_sc"
#> [4] "biomass_spr_sd_sc" "density_cod_rec_sc" "density_cod_sc"
#> [7] "density_fle_rec_sc" "density_fle_sc" "density_saduria_rec_sc"
#> [10] "density_saduria_sc" "depth_rec_sc" "depth_sc"
#> [13] "oxy_rec_sc" "oxy_sc" "sigma_E"
#> [16] "sigma_O" "temp_rec_sc" "temp_sc"
# I want the random effects to be dark gray...
pal <- brewer.pal(n = 8, name = "Dark2")
pal2 <- c(pal[1:5], "black", pal[6:8])
# Sort the terms so that random effects are at the top...
mfull_est <- mfull_est %>%
mutate(term2 = ifelse(term %in% c("sigma_E", "sigma_O"), 2, 1))
ggplot(mfull_est, aes(reorder(term, term2), estimate, color = group, fill = group, shape = scale)) +
geom_hline(yintercept = 0, linetype = 2, color = "gray40", alpha = 0.5) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0) +
geom_point(size = 2.5) +
scale_color_manual(values = pal2, name = "") +
scale_shape_manual(values = c(15, 19, 21, 17), name = "") +
scale_fill_manual(values = rep("white", 9), name = "") +
labs(y = "Estimate", x = "Standardized coefficient") +
scale_x_discrete(breaks = levels(mfull_est$term),
labels = c(expression(Herring[rec]),
expression(Herring[sub]),
expression(Sprat[rec]),
expression(Sprat[sub]),
expression(Cod[rec]),
expression(Cod[haul]),
expression(Flounder[rec]),
expression(Flounder[haul]),
expression(Saduria~entomon[rec]),
expression(Saduria~entomon[haul]),
expression(Depth[rec]),
expression(Depth[haul]),
expression(Oxygen[rec]),
expression(Oxygen[haul]),
expression(sigma[E]),
expression(sigma[O]),
expression(Temp[rec]),
expression(Temp[haul])
)) +
coord_flip() +
theme_plot() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.margin = unit(c(0.2, 0.2, 0.4, 0.2), "cm")) +
guides(color = "none", fill = "none", shape = guide_legend(ncol = 4))
ggsave("figures/effect_size.pdf", width = 17, height = 17, units = "cm")
# Just a test to see the labels were alright
ggplot(mfull_est, aes(reorder(term, estimate), estimate)) +
ggplot(mfull_est, aes(term, estimate)) +
geom_hline(yintercept = 0, linetype = 2, color = "red", alpha = 0.5) +
geom_point(size = 2) +
labs(x = "", y = "Standardized coefficient") +
coord_flip()
# Now save this for comparison in sensitivity analysis script
write.csv(mfull_est, "output/mfull_est.csv")
pred_map <- predict(mfull, newdata = pred_grid)
pal <- wes_palette("Zissou1", 21, type = "continuous")
plot_map_raster +
geom_raster(data = filter(pred_map, depth < 120 & depth > 10), aes(x = X*1000, y = Y*1000, fill = exp(est))) +
scale_fill_gradientn(colours = rev(pal), name = "Le Cren's condition factor") +
facet_wrap(~ year, ncol = 5) +
theme_facet_map() +
geom_sf(size = 0.1)
ggsave("figures/supp/condition/all_years_condition_map_covars.png", width = 6.5, height = 8.5, dpi = 600)
# And a smaller map for selected years
lab_size <- 1.7
plot_map_raster +
geom_raster(data = filter(pred_map, year %in% c(1994, 2001, 2008, 2018) & depth < 120 & depth > 10),
aes(x = X*1000, y = Y*1000, fill = exp(est))) +
scale_fill_gradientn(colours = rev(pal), name = "Le Cren's condition factor") +
facet_wrap(~ year, ncol = 2) +
theme(legend.key.height = unit(0.1, "cm"),
legend.key.width = unit(0.5, "cm"),
legend.position = c(0.9, .036),
legend.direction = "horizontal",
legend.background = element_rect(fill = NA),
legend.text = element_text(size = 6),
legend.title = element_text(size = 7.4)) +
geom_sf(size = 0.1) +
annotate("text", label = "Sweden", x = xmin2 + 0.25*xrange, y = ymin2 + 0.8*yrange, color = "black", size = lab_size) +
annotate("text", label = "Denmark", x = xmin2 + 0.029*xrange, y = ymin2 + 0.43*yrange, color = "black", size = lab_size, angle = 75) +
annotate("text", label = "Germany", x = xmin2 + 0.07*xrange, y = ymin2 + 0.025*yrange, color = "black", size = lab_size) +
annotate("text", label = "Poland", x = xmin2 + 0.55*xrange, y = ymin2 + 0.1*yrange, color = "black", size = lab_size) +
annotate("text", label = "Russia", x = xmin2 + 0.95*xrange, y = ymin2 + 0.2*yrange, color = "black", size = lab_size) +
annotate("text", label = "Lithuania", x = xmin2 + 1*xrange, y = ymin2 + 0.47*yrange, color = "black", size = lab_size, angle = 75) +
annotate("text", label = "Latvia", x = xmin2 + 0.99*xrange, y = ymin2 + 0.68*yrange, color = "black", size = lab_size, angle = 75)
ggsave("figures/condition_map_covars.pdf", width = 17, height = 17, units = "cm")
# Plot spatiotemporal random effect
plot_map_raster +
geom_raster(data = pred_map, aes(x = X * 1000, y = Y * 1000, fill = epsilon_st)) +
scale_fill_gradient2() +
facet_wrap(~ year, ncol = 5) +
theme_facet_map() +
geom_sf(size = 0.1)
ggsave("figures/supp/condition/epsilon_st_map.png", width = 6.5, height = 8.5, dpi = 600)
# Plot spatial random effect
plot_map_raster +
geom_raster(data = filter(pred_map, year == 1999), aes(x = X * 1000, y = Y * 1000, fill = omega_s)) +
scale_fill_gradient2() +
geom_sf(size = 0.1)
ggsave("figures/supp/condition/omega_s_map.png", width = 6.5, height = 6.5, dpi = 600)
# Fit a linear model to each prediction grid of the estimate over time
# https://community.rstudio.com/t/extract-slopes-by-group-broom-dplyr/2751/7
time_slopes_by_year <- pred_map %>%
mutate(id = paste(X, Y, sep = "_")) %>%
split(.$id) %>%
purrr::map(~lm(est ~ year, data = .x)) %>%
purrr::map_df(broom::tidy, .id = 'id') %>%
filter(term == 'year')
# Plot the slopes
time_slopes_by_year2 <- time_slopes_by_year %>%
separate(id, c("X", "Y"), sep = "_") %>%
mutate(X = as.numeric(X),
Y = as.numeric(Y))
plot_map_raster +
geom_raster(data = time_slopes_by_year2, aes(x = X * 1000, y = Y * 1000, fill = estimate)) +
scale_fill_gradient2(midpoint = 0, name = "Slope (condition factor~year)") +
geom_sf(size = 0.1)
<img src=“condition_model_cf_files/figure-html/calculate”spatial trend”-1.png” width=“672” style=“display: block; margin: auto;” />
ggsave("figures/supp/condition/spatial_trend_condition.png", width = 6.5, height = 6.5, dpi = 600)
# Prepare prediction data frame
nd_dep <- data.frame(depth_sc = seq(min(d$depth_sc), max(d$depth_sc), length.out = 100))
nd_dep <- nd_dep %>%
mutate(year = 2003L,
year_f = as.factor(2003),
ln_length_cm = 0,
biomass_her_sc = 0,
biomass_her_sd_sc = 0,
biomass_spr_sc = 0,
biomass_spr_sd_sc = 0,
density_cod_sc = 0,
density_cod_rec_sc = 0,
density_fle_sc = 0,
density_fle_rec_sc = 0,
density_saduria_sc = 0,
density_saduria_rec_sc = 0,
depth_rec_sc = 0,
oxy_sc = 0,
oxy_rec_sc = 0,
temp_sc = 0,
temp_rec_sc = 0
)
# Predict from full model
p_cond_dep <- predict(mfull, newdata = nd_dep, se_fit = TRUE, re_form = NA)
cond_dep <- ggplot(p_cond_dep, aes(depth_sc, exp(est),
ymin = exp(est - 1.96 * est_se),
ymax = exp(est + 1.96 * est_se))) +
geom_ribbon(alpha = 0.4) +
geom_line() +
theme(aspect.ratio = 1) +
labs(x = expression(Depth[haul]))
# Prepare prediction data frame
nd_oxy <- data.frame(oxy_sc = seq(min(d$oxy_sc), max(d$oxy_sc), length.out = 100))
nd_oxy <- nd_oxy %>%
mutate(year = 2003L,
year_f = as.factor(2003),
ln_length_cm = 0,
biomass_her_sc = 0,
biomass_her_sd_sc = 0,
biomass_spr_sc = 0,
biomass_spr_sd_sc = 0,
density_cod_sc = 0,
density_cod_rec_sc = 0,
density_fle_sc = 0,
density_fle_rec_sc = 0,
density_saduria_sc = 0,
density_saduria_rec_sc = 0,
depth_sc = 0,
depth_rec_sc = 0,
#oxy_sc = 0,
oxy_rec_sc = 0,
temp_sc = 0,
temp_rec_sc = 0
)
# Predict from full model
p_cond_oxy <- predict(mfull, newdata = nd_oxy, se_fit = TRUE, re_form = NA)
cond_oxy <- ggplot(p_cond_oxy, aes(oxy_sc, exp(est),
ymin = exp(est - 1.96 * est_se),
ymax = exp(est + 1.96 * est_se)))+
geom_ribbon(alpha = 0.4) +
geom_line() +
theme(aspect.ratio = 1) +
labs(x = expression(Oxygen[haul]))
# Prepare prediction data frame
nd_tem <- data.frame(temp_sc = seq(min(d$temp_sc), max(d$temp_sc), length.out = 100))
nd_tem <- nd_tem %>%
mutate(year = 2003L,
year_f = as.factor(2003),
ln_length_cm = 0,
biomass_her_sc = 0,
biomass_her_sd_sc = 0,
biomass_spr_sc = 0,
biomass_spr_sd_sc = 0,
density_cod_sc = 0,
density_cod_rec_sc = 0,
density_fle_sc = 0,
density_fle_rec_sc = 0,
density_saduria_sc = 0,
density_saduria_rec_sc = 0,
depth_sc = 0,
depth_rec_sc = 0,
oxy_sc = 0,
oxy_rec_sc = 0,
#temp_sc = 0,
temp_rec_sc = 0
)
# Predict from full model
p_cond_tem <- predict(mfull, newdata = nd_tem, se_fit = TRUE, re_form = NA)
cond_temp <- ggplot(p_cond_tem, aes(temp_sc, exp(est),
ymin = exp(est - 1.96 * est_se),
ymax = exp(est + 1.96 * est_se))) +
geom_ribbon(alpha = 0.4) +
geom_line() +
theme(aspect.ratio = 1) +
labs(x = expression(Temperature[haul]))
# Prepare prediction data frame
nd_cod <- data.frame(density_cod_sc = seq(min(d$density_cod_sc), max(d$density_cod_sc), length.out = 100))
nd_cod <- nd_cod %>%
mutate(year = 2003L,
year_f = as.factor(2003),
ln_length_cm = 0,
biomass_her_sc = 0,
biomass_her_sd_sc = 0,
biomass_spr_sc = 0,
biomass_spr_sd_sc = 0,
#density_cod_sc = 0,
density_cod_rec_sc = 0,
density_fle_sc = 0,
density_fle_rec_sc = 0,
density_saduria_sc = 0,
density_saduria_rec_sc = 0,
depth_sc = 0,
depth_rec_sc = 0,
oxy_sc = 0,
oxy_rec_sc = 0,
temp_sc = 0,
temp_rec_sc = 0
)
# Predict from full model
p_cond_cod <- predict(mfull, newdata = nd_cod, se_fit = TRUE, re_form = NA)
cond_cod <- ggplot(p_cond_cod, aes(density_cod_sc, exp(est),
ymin = exp(est - 1.96 * est_se),
ymax = exp(est + 1.96 * est_se))) +
geom_ribbon(alpha = 0.4) +
geom_line() +
theme(aspect.ratio = 1) +
labs(x = expression(Cod[haul]))
# Prepare prediction data frame
nd_spr <- data.frame(biomass_spr_sd_sc = seq(min(d$biomass_spr_sd_sc), max(d$biomass_spr_sd_sc), length.out = 100))
nd_spr <- nd_spr %>%
mutate(year = 2003L,
year_f = as.factor(2003),
ln_length_cm = 0,
biomass_her_sc = 0,
biomass_her_sd_sc = 0,
biomass_spr_sc = 0,
#biomass_spr_sd_sc = 0,
density_cod_sc = 0,
density_cod_rec_sc = 0,
density_fle_sc = 0,
density_fle_rec_sc = 0,
density_saduria_sc = 0,
density_saduria_rec_sc = 0,
depth_sc = 0,
depth_rec_sc = 0,
oxy_sc = 0,
oxy_rec_sc = 0,
temp_sc = 0,
temp_rec_sc = 0
)
# Predict from full model
p_cond_spr <- predict(mfull, newdata = nd_spr, se_fit = TRUE, re_form = NA)
cond_spr <- ggplot(p_cond_spr, aes(biomass_spr_sd_sc, exp(est),
ymin = exp(est - 1.96 * est_se),
ymax = exp(est + 1.96 * est_se))) +
geom_ribbon(alpha = 0.4) +
geom_line() +
#coord_cartesian(ylim = c(-4.65, -4.4)) +
theme(aspect.ratio = 1) +
labs(x = expression(Sprat[sub]))
# Prepare prediction data frame
nd_sad <- data.frame(density_saduria_rec_sc = seq(min(d$density_saduria_rec_sc),
max(d$density_saduria_rec_sc), length.out = 100))
nd_sad <- nd_sad %>%
mutate(year = 2003L,
year_f = as.factor(2003),
ln_length_cm = 0,
biomass_her_sc = 0,
biomass_her_sd_sc = 0,
biomass_spr_sc = 0,
biomass_spr_sd_sc = 0,
density_cod_sc = 0,
density_cod_rec_sc = 0,
density_fle_sc = 0,
density_fle_rec_sc = 0,
density_saduria_sc = 0,
#density_saduria_rec_sc = 0,
depth_sc = 0,
depth_rec_sc = 0,
oxy_sc = 0,
oxy_rec_sc = 0,
temp_sc = 0,
temp_rec_sc = 0
)
# Predict from full model
p_cond_sad <- predict(mfull, newdata = nd_sad, se_fit = TRUE, re_form = NA)
cond_sad <- ggplot(p_cond_sad, aes(density_saduria_rec_sc, exp(est),
ymin = exp(est - 1.96 * est_se),
ymax = exp(est + 1.96 * est_se))) +
geom_ribbon(alpha = 0.4) +
geom_line() +
#coord_cartesian(ylim = c(-4.75, -4.4)) +
theme(aspect.ratio = 1) +
labs(x = expression(Saduria~entomon[rec]))
# cond_dep + cond_oxy + cond_temp + cond_cod + cond_spr + cond_sad
# plot_annotation(tag_levels = 'A')
#
# ggsave("figures/supp/condition/conditional_effects.png", width = 6.5, height = 6.5, dpi = 600)
p_cond_dep2 <- p_cond_dep %>%
mutate(variable = "Depth (haul)") %>%
dplyr::select(est, est_se, depth_sc, variable) %>%
rename(value = depth_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)
p_cond_oxy2 <- p_cond_oxy %>%
mutate(variable = "Oxygen (haul)") %>%
dplyr::select(est, est_se, oxy_sc, variable) %>%
rename(value = oxy_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)
p_cond_tem2 <- p_cond_tem %>%
mutate(variable = "Temperature (haul)") %>%
dplyr::select(est, est_se, temp_sc, variable) %>%
rename(value = temp_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)
p_cond_cod2 <- p_cond_cod %>%
mutate(variable = "Cod (haul)") %>%
dplyr::select(est, est_se, density_cod_sc, variable) %>%
rename(value = density_cod_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)
p_cond_spr2 <- p_cond_spr %>%
mutate(variable = "Sprat (sub)") %>%
dplyr::select(est, est_se, biomass_spr_sd_sc, variable) %>%
rename(value = biomass_spr_sd_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)
p_cond_sad2 <- p_cond_sad %>%
mutate(variable = "S. entomon (rec)") %>%
dplyr::select(est, est_se, density_saduria_rec_sc, variable) %>%
rename(value = density_saduria_rec_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)
conds <- bind_rows(p_cond_dep2, p_cond_oxy2, p_cond_tem2,
p_cond_cod2, p_cond_spr2, p_cond_sad2)
ggplot(conds, aes(value, exp(est), ymin = exp(est - 1.96 * est_se),
ymax = exp(est + 1.96 * est_se))) +
geom_ribbon(alpha = 0.4) +
geom_line() +
facet_wrap(~variable, scales = "free_x") +
labs(y = "Le Cren's condition factor",
x = "Scaled value") +
theme_plot()
ggsave("figures/supp/condition/conditional_effects.png", width = 6.5, height = 6.5, dpi = 600)